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Abstract 



We show that the method of factorizing the evolution operator to fourth 
order with purely positive coefficients, in conjunction with Suzuki's method 
of implementing time-ordering of operators, produces a new class of powerful 
algorithms for solving the Schrodinger equation with time-dependent poten- 
tials. When applied to the Walker-Preston model of a diatomic molecule in a 
strong laser field, these algorithms can have fourth order error coefficients that 
are three orders of magnitude smaller than the Forest-Ruth algorithm using 
the same number of Fast Fourier Transforms. When compared to the second 
order split-operator method, some of these algorithms can achieve comparable 
convergent accuracy at step sizes 50 times as large. Morever, we show that 
these algorithms belong to a one-parameter family of algorithms, and that 
the parameter can be further optimized for specific applications. 

PACS: 31.15.-p, 02.70.Hm, 03.65.-W 
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I. INTRODUCTION 



Recently, we derived a new class of fourth order algorithms for solving both classical 
imj^ and quantum dynamical problems. These algorithms are based on factorizing the 
evolution operator e^*-"^"*"^-* to fourth order with purely positive coefficients and require know- 
ing the gradient of the force in the classical case and the gradient of the potential in the 
quantum case. The resulting algorithms are symplectic or unitary, respectively. While posi- 
tive coefficients are absolutely necessary for simulating the diffusion process in Monte Carlo 
algorithms @,|^, or doing imaginary time projections 0, they are not essential in quantum 
or classical algorithms. Nevertheless, we have shown that this class of gradient symplec- 
tic algorithms is far superior to existing fourth order algorithms with negative coefficients 
In this work, using Suzuki's method of implementing operator time-ordering, we 
prodouce a class of even more effective algorithms for solving quantum dynamical problems 
with explicit time-dependent potentials. Despite the vast literature on this subject |]§-|TB[, 
we believe our work has initiated a new direction in algorithm development. In the past, one 
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labors assiduously to avoid higher order commutators. Here, we show that their inclusion 
can yield algorithms of great efficiency. Our algorithm 4A, to be discuss below, is the fastest 
fourth order algorithm known, needing only four Fast Fourier Transforms (FFT) per itera- 
tion. This is only twice the computational effort of the second order split-operator method 
P JT^JT8[] , but the algorithm can converge at time step sizes an order of magnitude larger. 
Our optimized algorithms, when compared on an equal effort basis, have fourth order error 
coefficients that are three orders of magnitude smaller than Forest-Ruth's algorithm |T^ and 



a factor of 30 smaller than McLachlan's algorithm both are fourth order algorithms 
with negative coefficients. 

While these gradient symplectic algorithms are very efficient when the potential gradient 
is known analytically, they remain equally effective when the gradient is obtained numerically 
0. These algorithms are of particular interest in solving the time-dependent Schrodinger 
in a large 3D mesh. In 3D, even for a modest grid size of (256)^, the number of mesh points 
already exceed 16 millions. If the wave function array is double precision and complex, 
its storage alone would have required 268MB. For such a large number of grid points, any 
vector-matrix multiplication would be prohibitively expensive and must be avoided. For our 
algorithms, the costliest computational step is just the use of FFT. The unitary character 
and the large time step acceptance of these algorithms make them ideal for doing long time 
quantum simulations. 

The key problem in solving the Schrodinger equation with time-dependent potentials is 
the time ordering of operators. This problem is solved variously in the literature by trans- 
forming it into a classical problem [[l4| , p!5| , treating time as another "spatial coordinate" 
I^ , p!3|Jl6| , introducing auxiliary variables |T^, etc.. Most end up with some sort of time 



derivative operator, but none has the simplicity of Suzuki's method of directly implement- 
ing time-ordering via a forward time derivative operator. No time integration is necessary. 
Since this work is less accessible, but is of special relevance to the operator factorization 
approach of deriving algorithms, we summarize it in some detail in the next section. In 
Section III, we apply Suzuki's method and derive four gradient symplectic algorithms for 
solving the time- dependent Schrodinger equation. In Section IV, we use these algorithms to 



solve the Walker-Preston model |2T| and compare their convergent properties with existing 
algorithms. In Section V, we derive one-parameter families of these algorithms and show 
that they can be further optimized for specific applications. Section VI summarizes our 
conclusions. 



II. OPERATOR DECOMPOSITION OF ORDERED EXPONENTIALS 



For H{t) a time-dependent operator, the evolution equation 

d 



g^m = H{t)m, (1) 



has the operator solution 

rt+At 



il:{t + At) =T(exp H{s)ds)ij{t). (2) 

The time-ordered exponential not only has the conventional expansion 
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/•t+At i-t+At i-t+At i-si 

r(exp H{s)ds) = 1 + H{si)dsi + dsi ds2H{si)H{s2) + ■ ■ ■ , (3) 

but also the more intuitive interpretation 

rt+At 



T(exp / H{s)ds) = jim T(e^i:ti^(*+^^)), 

= lim e^^(*+^*) ■ ■ ■ e^^(*+^)e^^(*+^) . (4) 



There are many ways of solving the time-ordering problem. For this work on operator fac- 
torization, we prefer Suzuki's method 0, which directly implements time ordering without 
any additional formalism |]T2[ or auxiliary variables . 



Let D denotes the forward time derivative operator 

such that for any two time-dependent functions F{t) and G{t), 

F(t)e^*^G'(t) = F{t + At)G{t). (6) 

Suzuki's proved ^ that 

. ft+At . 

T(exp y H{s)ds) = exp[At{H{t) + D)]. (7) 



Using the more intuitive definition of time-ordering and invoking Trotter's formula, one 
proof only requires two lines: 

exp[At{H{t) + D)] = Jim (e^^We^^)", 

= lim e^^(*+^*) ■ ■ ■ e^^(*+^)e^^(*+^) (8) 

where property (|^) has been applied repeatedly and accumulatively. 

For the widely applicable case of H{t) = T + V{t), where only one of the operator is 
explicitly dependent on time, the short time evolution of can be written using ([^ as 

7/>(t + At)=e^*[^+^WVW, (9) 

which is just like the time-independent case but with an effective T = T + D. This suggests 
a two-step approach of deriving time-dependent algorithms. First, decompose e^*'"^"'"^''*^' in 
terms of e^*^ and e'^*^'-*^ using any factorization scheme applicable in the time-independent 
case. Next, since [T, D] = 0, factorize exactly 

and incorporate all time- dependent requirements by applying (P). For example, a second 
order factorization of (M) gives. 
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7<2)^ giAtTgAty(t)giAtf ^ 

= Qh^tT^AtV{t+At/2)^lAtT^ ^^^^ 

which is the well known midpoint algorithm for time-dependent problems. The other second 
order factorization gives the alternative second order algorithm, 

j-(2)^^lAtVit)^At¥^^AtVit)^ 

^ giAty(t)gAtCgAtrgiAty(t)^ 

^ giAty(t+At)gAtrgiAty(t)_ ^^2) 

Thus, for H{t) = T + V{t), the effect of time-ordering is to increment the time- dependence 
of each potential operator V{t) by the sum of time steps of all the T operators to its right. 

For the Schrodinger equation with a time-dependent potential, the wave function is 
evolved forward in a short time /St by 

^l,{e) = e^[^+^Wl^(0), (13) 
where e = -iAi, f + T,D ^iD, 

^=-^E^ V{t)^V{x,,t). (14) 

We will work in atomic units such that the kinetic energy operator has this standard form. 
Moreover, to do away with messy notations involving —i, we will use e as the time step 
variable everywhere. When e appears as the argument of the wave function or potential, it 
is to be understood that it denotes only the real time step variable At without the factor 
—i. (In this way, algorithms can be directly applied to the classical case with e purely real 
without change in form.) For conciseness, we will always regard the present time step as 
time zero. Thus, the two second order algorithms for solving the Schrodinger equation with 
step size At can be denoted simply as 

rf)(e)=e^^^e^^(^/2)e^^'^, (15) 
ri')(e)=e^^^(^)e^^e^^^(°). (16) 

Since the kinetic energy operator is diagonal in momentum space, the operator e^^ can 
be implemented as a vector- vector multiplication in Fourier space. Every occurrence of e^"^ 
requires two FFTs, one direct to Fourier space for the kinetic energy multiplication, and one 
inverse back to real space for the potential energy multiplication. To minimize the call for 
FFTs, one favors algorithms with the fewest occurrence of the kinetic energy operator. 

III. GRADIENT SYMPLECTIC FOURTH ORDER ALGORITHMS 

Following our two-step approach, we can transcribe any time-independent factorization 
algorithm into a time-dependent algorithm. For example, the well known Forest-Ruth (FR) 
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algorithm |T^(also discovered independently by Campostrini and Rossi [^,and Candy and 



Rozmus |2^) can now be transcribed to solve time- dependent problems as 
where, s = 2 -"^Z^, 

1 1 1 s , ^ 

-0 = vs = -j^^, .i = .2 = -2^, ^1=^3 = ^, ^2 = -^, (18) 

and Oi = ti, 02 = ti+t2- For easy identification, we adopt the convention of labelling the time 
step coefficients of operators T and V by and Vi respectively, and denote the intermediate 
time arguments of V by coefficients Oj. The coefficient of the first operator on the right will 
be denote by Vq or to, followed by paired coefficients of Vi and for i = 1,2, . . .n. If one 
were to decompose e'^*^^^^-' only in terms of e^^ and e^^, then Forest-Ruth is the only fourth 
order algorithm possible with 6 FFTs per iteration. The alternative algorithm with T and 
V interchanged, with appropriate modifications of the Oj coefficients, is also possible, but 
would have required 8 FFTs. Notice that some of the coefficients are negative, requiring 
backward propagation and evaluating the potential at a time prior to the present. This is 
a consequence of Suzuki's "no-go" theorem |^^, which proved that beyond second order, 
ge{T+v) Qannot be decomposed into a finite products of e*'^^ and e"^^^ with purely positive 
coefficients ti and Vi. Thus without exception, all higher order factorization algorithms 
heretofore proposed in the literature p|- pA| , p!5| -[T^ contain negative coefficients. 



It is the search for positive coefficient factorizations schemes that led one of us to 
derive fourth order symplectic algorithms for solving classical |I| and subsequently quan- 
tum dynamical problems. To circumvent Suzuki's "no-go" theorem, these factorization 
schemes employ an additional operator. 
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IVAT.Vll ^ l^i^J , (19) 

which is just the square of the gradient of the potential [^]. This can be computed an- 
alytically or numerically. For brevity, we will transcribe four previously derived gradient 
algorithms to their time-dependent form. In the next section, we will describe in more 
detail the one-parameter family of these algorithms. 

Our algorithm 4A (See also Ref. |]2B|), when applied to the time-dependent case of 
(IH), gives 

TW(e) = eiey(.)gi.Tg|ey(./2)gi.Tgiey(o)^ (20) 

with V defined by 

V{t)=V{t) + ^e'[V{t),[f,V{t)]], 

= V{t) + ^e'[V{t),[T,Vm. (21) 

Note that this is a crucial simplification, because [^(t), [-D, \^(t)]] = 0! Thus the addition 
of the forward time derivative operator D causes no additional complication to the gradient 
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term. ( This is very fortunate, because if one were required to keep the alternative double 
commutator [T, [l^(t),T]], this commutator would have given rise to complicated new terms 
involving dV{t)/dt and d'^V{t)/dt^ ! ) With the introduction of the gradient term, algorithm 
4A can achieve fourth order acccuracy with only four FFTs. Note also that in order to 
minimize the evaluation of the gradient term, the double commutator has been placed at 
the center. This choice seems obvious, but there is intrinsic freedom to redistribute the 
commutator term among the three potential operator without affecting the fourth order 
convergence of the algorithm. This is true for all algroithms described below. We will 
assume that this redistribution can be done when necessary. 
Similarly algorithm 4B can be transcribed as 

where 

4„ = <, = i(l--L), ;, = _L, a, = i(l--L). = + (23) 

and with V given by 

V{t) = V{t) + ^(2 - V?>)e'[V{t\ [T, Vm- (24) 
The time- dependent forms of algorithm 4C and 4D are respectively, 

r^'^\e) ^QbT^leV(5e/Q)^\eT^\eV{e/2)^\eT^leV(e/&)^\eT^ (25) 

r^4) ^ ^i,y(,)gi,rg|.v'(2./3)gi.Tg|ey(e/3)gi.rgi.y(o) ^ (26) 

where V{t) is as defined by (pT]). The number of FFT required by each algorithm is given 
in Table ffl. 



IV. SOLVING THE WALKER-PRESTON MODEL 

To demonstrate the effectiveness of these new algorithms, we use them to solve the Walker 
and Preston model ||2l| of a diatomic molecule in a strong laser field. Since this problem has 



been used by many authors |[T3|JT^JT7|JT8[| to test their time-dependent algorithms, it is an 
excellent choice for comparing our algorithms. The model is defined by the one dimensional 
Hamiltonian (in atomic units). 



H 



with 



1 92 



2/i dx^ 



+ V{x,t), 



(27) 



V{x,t) = 1/0(1 -e-"^)2 + Ax cos(cjt). 



(28) 
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and where /i = 1745, Vq = 0.2251, a = 1.1741, A = 0.011025, and uj = 0.01787. The wave 
function is initially chosen to be the Morse oscillator ground state, which has the corre- 



sponding ground state energy Eq = (a;o/2)[l — tuo/SVo] with uiq = ay2Vo//U. In conformity 
with the above authors, we discretize the wave function ilj{x,t) using 64 uniform points at 
Xk = —0.8 + kAx, with Ax = 0.08 . The natural time scale defined by the laser oscillation 
frequence is r = 27i/uj = 351.6 . 

In Fig.|l], the percentage errors of the total energy, 100|£'(1000r) — Econ\/ Econ, is plotted 
as a function of the step size At used in the calculation. Econ = 5.029155i?o is the converged 
value at very small time steps. The plotting symbols indicate calculated results. The lines 
are monomial fits in either (At)^ or (At)^. SO is our calculation using the second order split- 



operator algorithm ([Iq) . RS3 is Gray and Verosky's best convergent results [|I4[. They use 
a third order algorithm, but the error is degraded by the Magnus approximation to second 
order. Both can be well-fitted by a single quadratic 6j(At)^, with coefficients bso = 1-1 and 
bnss = 0.67. These are plotted as lines running through the data. The convergent range 
for both is below At < r/100 = 3.52. Algorithm RS3 requires three times as much effort 
as SO. For the same amount of computational effort, one can run SO three times at At/3 
and gain an reduction of (1/3)^ in its error coefficient. This projected convergence curve 
(650/9) (At)2 is plotted broken line and labelled as SO'. For the same effort (6 FFTs) 
one can also run the Forest-Ruth algorithm and obtain results shown as solid triangles, with 
a fitted line 9.7 x 10~^(At)^. This equal effort comparison clearly demonstrates the greater 
efficiency of fourth order algorithms. For errors in the range of 0.1 to 0.01 percent, FR's 
time step size can be 4 to 6 times as large as SO"s. This ratio would increase if greater 
accuracy is required. Algorithm 4A, which requires only 2/3 the effort of FR, can be fitted 
by 4.8 X 10~^(At)'^. This fit is plotted as a dotted line barely visible above the zero error 
line over the range of the plot. 

To discuss the convergence of our gradient symplectic algorithms, we greatly expanded 
the plotting range in Fig.^j. Here, we plot directly the convergence of E{1000t)/Eo. The 
plotting symbols indicate calculated results. We retained SO and FR for comparison. The 
SO result is well-fitted by Econ/ Eo + 0.054{At)'^ . Algorithm 4C and 4D yielded indistinguish- 
able results, despite the fact that 4D uses only 6 FFTs, two less than 4C. Obviously one 
should not use algorithm 4C in the present case; we only included it here for completeness. 
Since the FR algorithm is known to have rather large errors, we have also implemented 
McLachlan's fourth order algorithm EO] and obtained results labelled as M. This algorithm 



requires 8 FFTs per iteration and is the best algorithm tested by Sanz-Serna and Portillo [17 



Both FR and M are examples of fourth order algorithms with negative coefficients. The step 
size convergence of all fourth order algorithms can be very well-fitted by Econ/ Eq + di{AtY. 
These are plotted as lines going through data points. The coefficient \di\ for each algorithm 
is listed in Table |. Since the computational effort per time step is proportional to the 
number of FFTs, A^j, the error per unit effort, taking into account the variation of step 
size with effort, would be proportinal to \di\Nf, i.e., this is a measure of error of equal 
computational effort for all fourth order algorithms. Again, for example, algorithm 4A only 
requires half the number of of FFTs as McLachlan's algorithm. At a given At in running 
McLachlan's algorithm, one can execute algorithm 4A twice at At/ 2 and reduce its error 
by a factor of (1/2)'^. Thus despite the appearance in Fig.|^, algorithm 4A actually has a 
much smaller error per unit effort than McLachlan's algorithm. We normalize this equal 
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effort error to the value of FR and define the normahzed, equal-effort error coefficient as 
^eq = \di/dFR\{Ni/Npji)'^. This value for each algorithm is listed in Table |. Thus, excluding 
4C, our gradient algorithms are roughly a factor of 10^ smaller in error than FR and a factor 
of 3 to 8 smaller than McLachlan's algorithm. A even more useful measure is to consider 
Teff = ^eq^'^i which would give the time step size relative to RF for achieving the same error 
with equal effort. This is also given in Table |. From examining Fig.|I] and 0, it certainly 
seems reasonable that our algorithms can converge at time steps nearly 5 times as large 
as RF's and on the order of 30 times as large as SO"s. We will discuss further optimized 
algorithms 4ACB and 4BDA in the next section. 



V. ONE-PARAMETER FAMILY OF ALGORITHMS 

Algorithms 4A and 4C are special cases of the more general algorithm 

T^cb{(^) = e*=*^^e''^^^^''^^^e*2^'^e''2^^^"^^^e*^^'^e''^^^^"^''^e*°^'^ (29) 

By use of a Mathematica program [^,|^ to symbolically combine exponentials of operators, 
it is easy to determine these positive coefficients to yield a fourth order algorithm. For the 
above form of the algorithm, there is one free parameter, which we will take to be to- Given 
to, the rest of the coefficients are: 

ti=t2 = ]:-to, h=to, Vi = V3 = l- — t;2 = 1 - (fi +t^3), (30) 
/ o i — 2to 



where V here is given by 



V{t) = V{t) + ^e'[V{t),[T,Vm, (31) 

^2 



with ai = to, a2 = 1/2, a^, = 1 — to, and 

1 

Un = — 

12 



1 



1 - 2to 6(1 - 2to)=^ 



(32) 



For to = 0, one has algorithm 4A. For to = 1/6 ~ 0.167, one recovers algorithm 4C. Thus one 
can change continuously from algorithm 4A to 4C, and beyond, by "dialing" the parameter 
< to < |(1 — ~ 0.21. At the upper limit of to = |(1 — -^), V2 = 0, and the algorithm 
becomes a variant of algoirthm 4B, where the commutator term remained at the center 
rather than equally distributed between the two potential operators at positions vi and v^. 
We shall refer to this algorithm as 4B'. Algorithm 4B' can be continuously transformed to 4B 
by redistributing the commutator term from the center to both sides. For example, one can 
multiply the central commutator term by a factor 1 — a and add a/2 times the commutator 
term to each potential operator on the side. As a ranges from to 1, algorithm 4B' is 
continuously transformed to 4B. When a reaches 1, the central commutator disappears and 
the number of FFTs collapses from 8 to 6. Thus both algorithms 4A and 4B are singular 
end points of this general algorithm with discontinuous changes in the numbers of FFT. 
Algorithm 4B and 4D are special cases of the general algorithm 



8 



all requiring 6 FFTs. Here, we choose ti as the free parameter. The coefficients are then 

. . 6ti(ti-l) + l 1 

h = ti, t2 = l-2ti, V3 = Vo = V2=Vi = --Vo. (34) 



Here V is given by 



V{t) = V{t) + ^e'[V{t),[T,V{t)]], (35) 

V2 



with ai = ti, a2 = I — ti, and 



1 

" 48 



(36) 



Algorithm 4B corresponds to setting vo = and selecting the smaller of the two quadratic 
solution, ti = i(l — ^). When ti = 1/3, one obtains a variant of algorithm 4D, which we 
will denote as 4D'. Again, one can transform 4D' to 4D continuously by distributing the 
commutator term at positions Vi and V2 to Vq and V3. However, in order to keep the number 
of gradient terms at a minimum, we will not bother with this refinement here. For positive 
coefficients, we must have ^(1 — < ti < |. At the upper limit of ti = \, = 0, and 
the algorithm collapses back to algorithm 4A. Thus, there are two continuous families of 
algorithms, the C-type requiring eight FFTs and the D-type requiring six. They are joined 
at both ends by algorithms 4A and 4B with discontinuous changes in the numbers of FFTs. 
(Recently, Omelyan, Mryglod and Folk |28| have considered the one-parameter factorization 



schemes from 4A to 4C and 4B to 4D separately, without realizing that they can be further 
extended and joined into one.) 

Fig.0 shows that the fourth order convergence error is negative for algorithm 4A and 4B, 
and positive for 4C and 4D. Since algorithm 4A can be continuously changed to 4C, and 
algorithm 4B be transformed to 4D, there must be parameter values such that this fourth 
order error can be made to vanish. This immediately suggests a simple strategy for further 
optimization: for each application^ one can ffist minimized the convergence error for a short 
time run with respect to to in Tacb ^1 '^bda- can then use the algorithm at that 
value and at a large time step, to do long time simulations. Since the error term is generally 
time dependent, there is no guarantee that when it is small for one time, it will remain small 
for all time. For the important case of periodic time-dependence, one can hope to reduce 
the bound within which the error can fluctuate. This seems to work for the present case. 

From Fig. ^ as the algorithm is changed continuously from 4A to 4C and back to 4B, 
one should observe two zero error crossings. In varying to in Tacbi only see one. This 
suggests that the second crossing is when 4B' ^4B, which we did not consider in this work. 
The one zero crossing we observed is at to = 0.144, precisely between 4A(to=0) and 4C 
(to=l/6). A short run (t = lOOr) to determine the optimal value for to is shown as solid 
symbols with solid fitting lines in Fig.^. For to = 0.142 and to = 0.146, the convergence 
errors can be well fitted by fourth order lines as shown. At the value of to = 0.144, the 
fourth order coefficient is an order of magnitude smaller. We deem it sufficient to determine 
the zero-crossing parameter value to three-digit accuracy. The resulting long run is shown 
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as 4ACB in Fig.^. As shown in Table J, 5eq in this case is reduce by an order of magnitude 
and Te// lengthened to 8.5. The convergence error curve for 4ACB remained very flat in 
going from t = lOOr to t = lOOOr. 

In the case of T^'^a^ algorithm 4D' is sufficiently dissimilar to 4D such that the conver- 
gence error remained negative as 4B is morphed to 4D' at ti = 1/3. However, as ti increases 
above 1/3, the convergence error moves up and turns positive. Since 7^^^ ends up at 4A, 
the error must cross zero again on its way down. The first crossing is at ti = 0.35 . This is 
shown in Fig.^ with hollow symbol connected by fitted, broken lines. The second crossing 
is near ti = 0.46, but we do not bother to show this similar result. The long-run results 
using ti = 0.35 in Tg^^ is shown on Fig.^ as 4BDA. As listed in Table |I|, its equal effort 
error and effective time step factor are virtually identical to those of 4ACB. In both cases, 
without any additional computational effort, we have reduced the error by a factor of five 
and lengthen the time step size by a factor 1.5 . Both algorithms can achieve 10~^ accuracy 
in total energy at a time step size 50 times as large as SO'. 

A key advantage of these factorization algorithms is that one can obtain analytically 
their energy error terms by use of Mathematica. For fourth order gradient algorithms, there 
are four non-vanishing error operators in the Hamiltonian, 

AE(^) = (ei [VTVTV] + 62 [TTVTV] + 63 [VTTTV] + 64 [TTTTV] ) , (37) 

where are coefficient functions depending on to or ti, and [AiJC] = [A, [B, [C]]], etc.. 
The zero-crossing values of to in Tj^cB ^^"^ ^1 '^bda^ closely matched to the crossing 
points of the first two error coefficients, ei = 62- For 4ACB and 4BDA, the predicted 
zero-crossing values from solving ei — 62 = numerically are to = 0.14233 and ti = 0.35023 
respectively, in excellent agreement with empirical values. It seems that out of the four error 
operators, only two are dominant and they are opposite in sign. It should be noted that 
these crossing points have nothing to do with the minimum of Yli=iS^- These algorithms 



are not optimized by simply minimizing the sum of squares of error coefficients . Further 
discussions on the above one-parameter family of gradient algorithms will be presented in a 
separate publication. 



VI. CONCLUSIONS 

In this work, we have derived a new class of fourth order algorithms for solving the 
time-dependent Schrodinger equation with explicit time-dependent potentials. This class of 
algorithms is characterized by having only positive factorization coefficients and can achieve 
great efficiency by knowing the gradient of the potential. For the Walker-Preston model, on 
an equal effort basis, the convergence errors of these algorithms can be 5000 and 30 times 
smaller than negative-coefficient algorithms such as the Forest-Ruth and the McLachlan 
algorithm, respectively. These gradient algorithms should be further tested on more realistic 
problems where the gradient of potential may have to be computed numerically. 

Our discovery of the one parameter family of gradient algorithms illustrates the power 
of the operator factorization approach of solving any evolution equation. The ability to tai- 
lor a specific algorithm for a particular application refiects great versatility of our method. 
Our work encourages further systematic study of these and other parametrized families of 



10 



algorithm and their optimization. In particular, one should explore the freedom in redis- 



tributing the commutator term |23] in going from 4B to 4B' and from 4D to 4D'. We have 
not done so here in order to concentrate on algorithms with the minimum number of gradi- 
ent evaluations. With the development of these powerful gradient algorithms for solving the 
Schrodinger equation directly, we can see no practical advantage in solving the correspond- 
ing classical problem using classical methods. This is in agreement with Sanz-Serna and 
Portillo's earlier assessment [IT? . 



This work spurs further interest in finding six and higher order factorization algorithms 
with purely positive coefficients. Despite intense effort, we have yet to find a sixth order 
factorization scheme with positive coefficients. It is most likely that one doesn't exist, even 



with the inclusion of the double commutator term. A recent work reached the same 
conclusion, but offered no proof either. Such a proof, if exists, would make these fourth 
order gradient algorithms unique. There would be no higher order generalizations. This 
would raise many other interesting questions, such as what other operators are necessary 
for a higher order positive time step factorization, what would be the optimal sixth order 
algorithm with minimal of negative coefficients, etc.. All such investigations would deepen 
our understanding of the operator factorization method of solving evolution equations. 
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FIGURES 




FIG. 1. The percentage energy error at t = lOOOr for various algorithms in solving the 
Walker-Preston model. The plotting symbols are results of calculations. The lines are second 
or fourth order fits. SO is the split-operator algorithm (|l^), SR3 is one of Gray and Verosky's 
algorithm, FR is the Forest-Ruth algorithm (17), and 4A is the gradient algorithm (pO|). SO' is an 
equal computational effort version of SO for comparing with SR3 and FR. See text for details. 
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FIG. 2. The long time energy convergence error for various algorithms in solving the 
Walker-Preston model. Eq is the Morse oscillator's ground state energy. M is McLachan's fourth 
order algorithm. The gradient algorithms 4A, 4B, 4C, 4D are defined by (l20|),(p^),(^5|), and (|26|) 
respectively. The optimized one-parameter algorithms 4ACB and 4BDA are described by (p9|) and 
(|3^ . See text for details. 
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FIG. 3. The short time convergence error as a function of the parameters to and ti. The sohd 
symbols fitted by sohd hnes are results for algorithm 4ACB at three values of to. The hollow 
symbols with broken lines are results for algorithm 4BDA at three values of ti . 
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TABLES 



TABLE I. Equal computational effort comparison of all fourth order algorithms discussed in 
this work. All algorithms except FR and M are our new time-dependent algorithms. Ni is the 
numbers of FFTs per execution for each algorithm, di is the corresponding fourth order error 
coefficient, Sgq is the equal-effort fourth order error coefficient normalized to FR's value, and Tg// 
is the effective time-step size relative to FR's value, e.g., for the same amount of effort, algorithm 
4ACB can achieve the same convergence error of FR at a time step 8.5 times as large. 





FR 


M 


4A 


4B 


4C 


4D 


4ACB 


4BDA 
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8 


4 


6 


8 


6 


8 


6 


\di\ 


5.0 X 10-5 


1.3 X IQ-'^ 


2.4 X lO"'^ 


0.5 X 10"^ 


1.0 X lO"'^ 


1.0 X lO"'^ 


3.0 X 10-9 


9.5 X 10-9 




1 


8.2 X 10-3 


0.95 X 10-3 


1.0 X 10-3 


6.3 X 10-3 


2.0 X 10-3 


1.9 X 10-^ 


1.9 X 10-^ 


^eff 


1 


3.3 


5.7 


5.6 


3.5 


4.7 


8.5 


8.5 
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